cls
clear all
set more off
capture log close 
version 16
log using Aggregation_Cases1234_19March_2022.log, replace
*This was created by Aggregation_Cases1234_19March_2022.do
qui local rep = 500
qui local n=100000

foreach sigmax of numlist  0.20  { 
foreach sigma  of numlist  0.40  { 
di
di "----------------------------------------------------------------------"
if `sigma'==0&`sigmax'==0 di "CASE 1"
if `sigma'!=0&`sigmax'==0 di "CASE 2"
if `sigma'==0&`sigmax'!=0 di "CASE 3"
if `sigma'!=0&`sigmax'!=0 di "CASE 4"
di "----------------------------------------------------------------------"
di
foreach dgp of numlist 2 { // 1 2 Just use case 2 so that PPML is not optimal
foreach t of numlist 60 { 

clear 
set seed 24022022
qui local obs=`n'*`t'
qui local smpl=max(`obs',`rep')
qui set obs `smpl'
qui g obsr=_n
qui local A=5 
qui local B=2

qui gen ID=1+int((_n-1)/`t') if _n<=`obs'
qui gen T=1+`t'*(((_n-1)/`t')-int((_n-1)/`t')) if _n<=`obs'
qui xtset ID T

qui g bp=. if _n<=`rep'
qui g bg=. if _n<=`rep'
qui g wb0=. if _n<=`rep'

forvalues r = 1/`rep' {
qui g beta=rnormal(0,`sigma') if _n<=`obs'
qui bysort ID: egen betax=mean(beta) if _n<=`obs'
qui replace beta=-1+betax*sqrt(`t')

qui g x=rnormal() if _n<=`obs'
qui bysort ID: egen mx=mean(x) if _n<=`obs'
qui replace x=mx*sqrt(`t')
drop mx
if `sigmax'!=0 qui replace x=rnormal(x,sqrt((1-`sigmax')*(x<0)+(1+`sigmax')*(x>0))) if _n<=`obs'
qui bysort ID: egen mx=mean(x) if _n<=`obs'

if `dgp'<1.5 qui g double y=rchi2(rnbinomial(exp(-ln(`t')+beta*x+0.4)/(`A'-1),1/`A')) if _n<=`obs'
if `dgp'>1.5 qui g double y=rchi2(rnbinomial(1/(`B'*`t'),1/(1+`B'*exp(beta*x+0.4))))  if _n<=`obs'
qui replace y=0 if y==.&_n<=`obs' 
qui bysort ID: egen my=sum(y) if _n<=`obs'
sort obsr

qui g double yb0=y*beta if _n<=`obs'
su yb0 if _n<=`obs', meanonly
local mu0=r(mean)
su y if _n<=`obs', meanonly
qui replace wb0=`mu0'/r(mean) in `r'

qui glm my mx if T==1&_n<=`obs', irls robust link(log) family(poisson)
qui replace bp= _b[mx] in `r'
mat b=e(b)
qui glm my mx if T==1&_n<=`obs', robust link(log) family(gamma) from(b)
if e(converged)==1 qui replace bg= _b[mx] in `r'

drop x y* mx my* beta*  
}
di
di "DGP = `dgp', S = " `sigma' ", Sx = " `sigmax' ", t = " `t' ", n = " `n' ", 
su bp bg wb0 if _n<=`rep'
di "----------------------------------------------------------------------"
di "----------------------------------------------------------------------"
di
}
}
}
}
capture log close
exit
